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Massive Black Hole recoil in high resolution hosts 
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The final inspiral and coalescence of a black hole binary can produce highly beamed gravitational wave radiation. To 
conserve linear momentum, the black hole remnant can recoil with "kick" velocity fkick < 4000 km/s. We present two 
sets of full N-body simulations of recoiling massive black holes (MBH) in high-resolution, non-axisymmetric potentials. 
The host to the first set of simulations is the main halo of the Via Lactea I simulation (Diemand et al. 2007). The nature 
of the resulting orbits is investigated through a numerical model where orbits are integrated assuming an evolving, triaxial 
NFW potential, and dynamical friction is calculated directly from the velocity dispersion along the major axes of the main 
halo of Via Lactea I. By comparing the triaxial case to a spherical model, we find that the wandering time spent by the 
MBH is significantly increased due to the asphericity of the halo. For kicks larger than 200 km/s, the remnant MBH does 
not return to the inner 200 pc within 1 Gyr, a timescale an order of magnitude larger than the upper limit of the estimated 
QSO lifetime. The second set of simulations is run using the outcome of a high-resolution gas-rich merger (Mayer et al. 
2007) as host potential. In this case, a recoil velocity of 500 km/s cannot remove the MBH from the nuclear region. 

© 2008 W1LEY-VCH Verlag GmbH & Co. KGaA, Weinheim 



1 Introduction 

In the context of the currently favored ACDM cosmogony, 
where large halos are assembled through the hierarchical 
merging and accretion of small progenitors, MBH mergers 
should be common. While dynamical friction against the 
dark matter and baryons can only form binaries with separa- 
tions of ~1 pc, other dynamical process such as stellar ejec- 
tion and gas drag can efficiently reduce the binary orbital 
separation down to ~0.001 pc, the regime where gravita- 
tional wave radiation dominates the orbital energy loss. This 
radiation is typically anisotropic due to asymmetries associ- 
ated with the black holes' mass and spin, causing the center 
of mass of the system to recoil in order to balance the linear 
momentum carried away by the gravitational wave radiation 
(Bekenstein 1973; Fitchett & Detweiler 1984; Favata et al. 
2004). The "kick" velocity of the remnant depends on the 
mass ratio (m^/mi < 1) and the spin parameters (ai, 0,2) 
of the binary, but not on the total mass of the system. 
Early estimates of the recoil velocity, framed in the post- 
Newtonian regime for non-spinning, unequal mass black 
holes (Fitchett 1983; Redmount & Rees 1989), yielded ve- 
locities in the range 100 < Ukick < 500 km/s which were 
successfully reproduced by numerical data (Baker et al. 
2006b; Gonzalez et al. 2007b; Herrmann et al. 2007). Simu- 
lations with varying (arbitrary) spin orientations (Campan- 
elli et al. 2007b) and mass ratios (Baker et al. 2008) show 
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that recoil velocities can be significantly larger, reaching 
Wkick ~ 2000 km/s (Campanelli et al. 2007a; Gonzalez et 
al. 2007b) and are predicted to be as large as Wkick ~ 4000 
km/s for unequal mass (ma/mi = 1/3), maximally spin- 
ning black holes (Campanelli et al. 2007a). In order to ad- 
dress the issue of whether off-nuclear QSO/AGN are com- 
mon and detectable, we ought to understand the dynam- 
ics of their orbits. The radial orbit of a recoiling MBH in 
a spherically symmetric potential was studied analytically 
by Madau & Quataert (2004). They showed that large kicks 
(400 km/s) can displace MBHs a few tens of kiloparsecs 
away from the center of a Milky-Way size stellar bulge. 
After the kick, the MBH undergoes several oscillations be- 
fore decaying back to the bottom of the potential. Most of 
the orbital energy is lost during the MBH passages through 
the center, where dynamical friction is most efficient. Gua- 
landris & Merritt (2008) substantiated these results by per- 
forming direct summation N-body simulations of MBH re- 
coil in a spherically symmetric galaxy, modeled as a binary- 
depleted core-Sersic profile (Graham et al. 2003). They find 
that after reaching the core radius, the MBH and the core ex- 
periences damped oscillations with bulge's center of mass, 
which decelerate the remnant MBH until it reaches thermal 
equilibrium with the surrounding stars. They estimate that 
this final stage may take up to 1 Gyr in large galaxies. As we 
show in Sect. 12.21 the MBH wandering time can be largely 
increased if the host potential is triaxial, while the presence 
of gas can significantly damp the overall MBH motion. 
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Fig. 1 (online colour at: www.an-journal.org) Resulting 
3D orbit corresponding to simulation VL080. The first 0.5 
Gyr are plotted in yellow, the following 0.5 Gyr in orange 
and the remaining 0. 1 Gyr is plotted in dark red. The box 
has side length L — 1.2 kpc. 

After a major merger, however, the remnant dark matter 
halo tends to be prolate (e.g. Novak et al. 2006), a factor that 
can play a major role for large kicks. 

In the following we carry out N-body simulations of re- 
coiling MBH in two high-resolution non-axisymmetric po- 
tentials: the Via Lactea I simulation's main halo, a triaxial 
dark-matter-only host with known triaxiality (Kuhlen et al. 
2007) and evolution parameters (Diemand et al. 2007), and 
the gas-rich galaxy merger remnant described in Mayer et 
al. (2007). In the dark-matter only case, we present a numer- 
ical model that successfully reproduces the main features 
of the MBH orbits, and can therefore be used to estimate 
the fallback time and maximum displacement for any given 
kick velocity (Guedes et al. 2008, in preparation). 

2 Recoil in a triaxial dark matter halo 
2.1 Simulations 

Simulations of recoil in a triaxial dark matter halo are run 
using the entire volume of the Via Lactea I simulation (here- 
after VLI), a periodic box of size L = 90 Mpc. The high- 
resolution region contains over 230 million particles of mass 
of m p = 2.1 x 10 4 M0. The exquisite resolution of VLI al- 
lows us to adopt the mass of SgrA* M. = 3.7 x 1O 6 M 
(Ghez et al. 2005) for the MBH and a force softening length 
e = 90 pc. The five halo + MBH systems are evolved for 1 . 1 
Gyr using PKDGRAV (Stadel 2001), a tree algorithm that 
includes up to hexa-decapole moments to reach high accu- 
racy in the force calculation. The five runs correspond to 
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Fig. 2 Simulations VL080-VL400. Resulting orbits of re- 
coiling MBHs from A^-body simulations in the VLI main 
halo. 

five kick velocities w kick = 80, 120, 200, 300, and 400 km/s 
for runs labeled VL080 to VL400 respectively. 

We place MBH at the densest point of the halo, with 
phase-space coordinate w c , at an initial redshift z\ = 1.54, 
300 Myr after the halo suffered its last major merger. At 
this epoch the halo is characterized by i?2oo = 187 kpc and 
M200 = 1.02 x 1O 12 M , where JS200 is defined as the ra- 
dius within which the enclosed average density is 200 times 
the critical density. The MBH particle is tracked at every 
time-step and its position and velocity are measured with 
respect to w c . 

A 3D rendition of the orbit for simulation VL080 is 
shown in Fig. Q] Orbits for simulations VL080-VL400 are 
shown in Fig. [2] These orbits sample a large volume of VLI 
and reveal the asphericity of the main halo which causes an- 
gular momentum transfer to the y and z components of the 
MBH velocity. As will be shown in more detail in Guedes 
et al. (2008, in prep.), this increases the MBH wandering 
time, especially in cases where the kick is large enough 
to leave the baryonic component behind. For simulation 
VL080, which features the smallest kick velocity fkick = 80 
km/s, we are able to obtain both a maximum displacement 
of 1 .2 kpc and a fall back time of 1 . 1 Gyr and therefore it can 
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be used to constrain the semi-analytical model described on 
Sect. EH 



2.2 Triaxial model 

The orbit of a MBH can be characterized by a conservative 
force due to the dark matter potential V$, and a damping 
dynamical friction term. This generates a system of six cou- 
pled differential equations that can be separated as follows: 

r = v, (1) 

v = -V$ + / DF . 

The triaxial dark matter potential is modeled as modified 
NFW profile (Navarro, Frenk & White 1997), 

GM 200 m(l+x T ) 
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where — ?'t /?"s> with r s the scale radius, and p and q are 
the triaxiality parameters. 

The classical Chandrasekhar dynamical friction formula 
is not valid in a triaxial system, because the velocity dis- 
persion is non-isotropic and the velocity distribution devi- 
ates from Maxwellian. We adopt the Pesce et al. (1992) 
generalization of the dynamical friction formula for non- 
axisymmetric systems, 

/df = -riviei - T 2 v 2 e 2 ~ T 3 v 3 e 3 , (4) 
where Vi is the component of the black hole velocity along 
the principal axis e; of the halo's velocity dispersion ten- 
sor, and Ti are the dynamical friction coefficients given by 
(Pesce et al. (1992);Vicari et al.(2007)) 
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with €i = (Tij <Ji and o\ is the largest eigenvalue. 

The six components of the symmetric velocity disper- 
sion tensor, defined as er?- = (viVj) — (vi)(vj) are mea- 
sured as a function of radius for every output of the VLI 
simulation. To get the three principal velocity dispersions, 
we diagonalize afj, obtaining the principal eigenvalues a\, 
<j|, and erf as a function of radius. To simplify our calcula- 
tions, we perform analytical fits to a\, 
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where the best fit parameters A, B,C, D,n,m vary as a 
function of redshift. 

Several studies have shown that a constant value for the 
Coulomb logarithm, In A, does not match simulations re- 
sults. We adopt the Maoz (1993) formalism 
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Fig. 3 (online colour at: www.an-journal.org) Compari- 
son of spherical and triaxial models of the MBH orbits after 
suffering a gravitational wave recoil. The N-body simula- 
tion results (red) are superposed to the numerical model re- 
sults (orange). 



where po is the interior density, i? max = r s , and we solve 
for r m in = GM,/a 2 (r), the sphere of influence of the 
MBH. The resulting orbits are shown in Fig. [3] (right). Con- 
trary to the spherical model (left), the orbits in a triaxial 
halo are non-radial, which keeps the MBH orbiting away 
from the center of the host, reducing the mean contribution 
of dynamical friction to the motion and extending the wan- 
dering time. This model qualitatively reproduces the results 
of the N-body simulations described in Sect. 12.11 

3 Recoiling MBHs in a gas merger 

We present the preliminary results of a simulation of a 
recoiling MBH of mass M. = 3.6 x 1O 7 M in a multi- 
component galaxy using the SPH + N-body code GASO- 
LINE (Wadsley, Stadel & Quinn 2004). The host is the 
merger remnant of a high-resolution simulation by Mayer et 
al. (2007), which features a gas particle mass of m p = 3000 
Mq, a gravitational softening e = 2 pc, and a central MBH 
in each of the progenitor galaxies. The nuclear region of 
the merger remnant contains 2 million particles, contains a 
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Fig. 4 (online colour at: www.an-journal.org) Projected 
density plot of the nuclear region of the merger remnant. 
A wandering MBH is seen at apocenter as an over-density 
near the nuclear region at t = 5.4 Myr after suffering a 500 
km/s recoil kick. The image corresponds to an area of 30 pc 
x30 pc. 

mass of M g — 3 x 10 9 Mq and a density that ranges from 
p = 10 _2 -10 5 atoms/cm 3 . By the end of their simulation, 
the MBHs have formed a binary with semi-major axis equal 
to the softening length of the simulation. 

We take the simulation from this point forward. We re- 
move both MBHs and evolve the simulation until the nu- 
cleus becomes smooth, place our MBH at the center of mass 
of the system, and kick it with recoil velocity ^kick = 500 
km/s. In this case, the MBH was displaced only 30 pc from 
the center of mass. A simple Bondi accretion estimate yields 
an accretion rate of M acc = 0.06 M Q /yr, and therefore the 
MBH can become an active AGN during its wandering time. 

As seen in Fig. [4] the MBH produces an over-density 
signature on the nuclear gas at the apocenter of its orbit, a 
feature that can be associated with the emission of X-ray 
radiation (e.g. Fujita 2008). 

4 Discussion 

If MBH mergers are common in the context of ACDM, then 
the "rocket effect" should also be common, especially at 
high-redshift where the bulk of the mass assembly occurs. 
Because recoiling MBHs carry gas and stars with them, the 
detection of off-nuclear QSOs is possible by finding veloc- 
ity shifts between the narrow line region associated with gas 
accreting onto the MBH, and the broad line emission associ- 
ated with the galaxy left behind (Bonning & Shields 2007). 
The best off-nuclear QSO candidate today is thought to be 
powered by a recoiling SMBH of mass M, = 6 x 10 s M 
which travels at a velocity of 2650 km/s with respect to the 
host galaxy (Komossa et al. 2008). 



In order to assess the detectability of so-called naked 
QSOs, we ought to understand the dynamics of the recoil- 
ing MBH orbits and the environment provided by the host 
potential. While detectability requires the presence of an ac- 
cretion disk around the remnant MBH, the medium of its 
host galaxy can damp its motion significantly, disfavoring 
its detection. However, if the kick is large enough that, while 
keeping an accretion disk and stellar sphere on influence, 
the MBH can reach as far as the (generally) prolate dark 
matter halo, dynamical friction is much less efficient and 
therefore the detection probability is mostly determined by 
the lifetime of the QSO. 
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